/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TAMOMDIS
  \brief This is a global class, to take control of the whole flow of the
  program. It is also responsible for generating various momentum distributions.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/08
  \date Last modified: 2022/04/19 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include <string>
#include <vector>
#include <complex>
#include <fstream>
#include <unistd.h>

using std::ofstream;
using std::string;
using std::vector;

typedef std::complex<double> cdouble;

#ifndef _TAMOMDIS_h_
#define _TAMOMDIS_h_

class TABound;
class TASMatrix;
class TAMOMDIS_M;
class TAConvolve;
class TANucleus;

class TAMOMDIS{
public:
  TAMOMDIS(int argc, char *argv[]); ///< the general constructor
  /// the constructor for single-configuration mode
  TAMOMDIS(const string &cfgfile, bool isconv = false,
    bool ismom = false, bool isbound = false, bool exCS = false);
  virtual ~TAMOMDIS();

  void GetOpt(int argc, char *argv[]); ///< parse command options
  void PromptHelp(); ///< show manual for input rules
  void Go(); ///< here we go

  bool IsProtonRemoval() const; ///< -1: undefined; 1: true (zv: 1); 0: false (zv: 0)
  void Configure(); ///< solve Rl, generate and fit S-matrix
  void SetExCS(bool opt = true){ fExCS = opt; }
  void Parallel(); ///< calculate dsigma/dkz
  double Total(); ///< calculate total knockout cross sections
  TABound *GetBound();
  double Rl(double r);
  double u(double r);
  TASMatrix *Sc(); ///< \retval the S-matrix object of the core-target
  TASMatrix *Sv(); ///< \retval the S-matrix object of the valence-target
  cdouble Sc(double b);
  cdouble Sv(double b);

  /// convolve the proj. c.m. momdis (m) to give the exp. distri. stored in fout
  /// \param cc: normalization constant
  /// \param fconv: convolved f, return value of Convolve()
  void Convolve(int l, double *x, double *f, double cc, vector<double> &xconv,
    vector<double> &fconv, ofstream &fout);

  // getters for the outputs //
  int nmom() const{ return fNmom; }
  int nmomconv() const;
  bool IsConv() const{ return fIsConv; }
  double kz(int i) const;
  double kzconv(int i) const;
  double str(int m) const{ testm(m); return fSigmaStr_M[m]; }
  double dif(int m) const{ testm(m); return fSigmaDiff_M[m]; }
  double tot(int m) const{ testm(m); return fSigmaTotal_M[m]; }
  double str() const{ return fSigmaStr; }
  double dif() const{ return fSigmaDiff; }
  double tot() const{ return fSigmaTotal; }
  const vector<double> &momstr(int m) const{ testm(m); return fMomStr_M[m]; }
  const vector<double> &momdif(int m) const{ testm(m); return fMomDiff_M[m]; }
  const vector<double> &momtot(int m) const{ testm(m); return fMomTotal_M[m]; }
  const vector<double> &momstr() const{ return fMomStr; }
  const vector<double> &momdif() const{ return fMomDiff; }
  const vector<double> &momtot() const{ return fMomTotal; }
  const vector<double> &momstrConv() const{ return fMomStrConv; }
  const vector<double> &momtotConv() const{ return fMomTotalConv; }
  string banner() const{ return fBanner; }

  friend class TAMOMDIS_M;

protected:
  void testm(int &m) const; // test the validity of m

  /// the component to calculate momentum distribution of the core
  TABound *fBound;
  TASMatrix *fSc, *fSv;
  int l; ///< the angular momentum quantum number
  int fNmom; ///< number of momentum bins
  TAMOMDIS_M *fMOM_M; // calculate m-specific c.s.
  bool fHasConfigured; // whether Configure() has been called or not

  /// c.s. components of different angular momentum projections (m)
  vector<double> fSigmaStr_M, fSigmaDiff_M, fSigmaTotal_M; ///< [l+1]
  double fSigmaStr, fSigmaDiff, fSigmaTotal;
  /// -m and m have the same dsigma/dkz. Their sum ([-m]+[m]) are stored for m!=0
  vector<vector<double>> fMomStr_M, fMomDiff_M; ///< [l+1][kNmom]
  vector<vector<double>> fMomTotal_M; ///< sum of str and diff momdis: [l+1][kNmom]
  vector<double> fMomStr, fMomDiff, fMomTotal; ///< [kNmom]
  /// convolved with exp resolutions: [kNmom]
  vector<double> fMomConv, fMomStrConv, fMomTotalConv;

  /// convolve the proj. c.m. momdis to give the exp. distri. //
  TAConvolve *fConv;
  /// true: only solve the bound state; false(default): going through all procedures
  bool fIsBound;
  bool fIsMom; ///< true: calculate momdis; false(default): calculate total c.s.
  bool fIsBatch; ///< true: batch mode; false(default): single configuration mode
  int fIsProtonRemoval; ///< -1: not defined; 1: true (zv == 1); 0: false (zv == 0)
  bool fExCS; ///< skip calcualtion of bound and Smatrices when true
  string fCfgFile; ///< the config file
  bool fIsConv; ///< true: do the convolution with exp resolutions; false(default): not
  vector<TANucleus *> fnv; ///< for batch mode, each element calculates a nulceus
  string fExpResCfgFile;  ///< yaml cfg file containing the exp resolutins
  string fBanner; ///< the reaction info
};

#endif
